用图层叠加方法绘制环形进化树
1. 开发背景
系统发育树可以应用于群体遗传学，生态学以及流行病学等研究领域。随着相应研究的深入，我们往往需要将进化树与相关数据联系起来。比如生态学中的，末端节点的物种分类，丰度等等信息。通过这样的联系也便于比较研究。目前，已有一些网页软件被开发，比如iTOL
(Letunic and Bork 2019), EvolView (Subramanian et al.
2019)。使用这些web工具仍然很难将外部数据与进化树进行整合。前几年的一款基于python开发的GraPhlan
(Asnicar et al.
2015)，可以实现在环形进化树外圈添加相关数据信息。但其支持的几何图层有限，而且一开始的配置文件对于大多数研究者都很难构建，其支持的树文件格式也较少，支持的可视化进化树的布局也只有circular。相比之下，本实验室开发的ggtree(Yu et al. 2017, 2018; Yu
2020)，是基于图形语法开发的R包，使用方便，支持多种树文件输入，有利于可重复计算，也很容易整合到现有分析流程中，自其发表以来受到相关研究领域的热门关注。目前ggtree可以通过geom_facet或者facet_plot来将phylogenetic tree与一些相关的metadata数据可视化联系起来。但其支持的进化树布局也较少。对于circular等layout还不能通过图形语法来简单进行图层叠加。为了填补这个空白，最近，我们开发了一个R包ggtreeExtra。它是基于图形语法，可以用ggplot2的mapping等简单操作将外部数据可视化成合适图层,然后叠加到进化树上，而且支持多种布局的进化树。
2. 安装与使用
ggtreeExtra的开发地址点文章左下角的阅读原文即可。这个包还未发布到CRAN或者Bioconductor上。所以安装可以通过github来。
if (!require("remotes", quietly=TRUE)){
install.packages("remotes")
}
remotes::install_github(YuLab-SMU/ggtreeExtra)
因为ggtreeExtra是通过图形语法来写的，所以使用方法与ggplot2类似。主要是通过geom_fruit来进行图层叠加，该函数的几个参数解释。其中data为绘制图层的数据，数据格式与你需要指定绘制图层的数据格式一致。geom参数则为几何图层指定。目前支持，geom_bar,geom_point, geom_tile, geom_star, geom_boxplot,geom_violin。mapping参数为aes映射参数与ggplot2一致。offset参数为控制不同图层的间距大小，其代表该距离与树的比例，，默认0.03，则为树的大小0.03倍。pwidth参数控制该图层的宽度，其代表该图层与树的比例，默认为0.2，即树为1，该图层为0.2。还有...参数，用来传递指定的几何图层的参数。
3. 案例
3.1 案例1
树形图，外圈叠加连续型热图，binary热图与柱状图。在本例子除了树的信息其他的外部数据均为随机生成。而且这些外部数据的格式均为ggplot2相应的几何图层的数据输入格式。
#图层叠加的包
library(ggtreeExtra)
#绘图
library(ggplot2)
#绘制进化树
library(ggtree)
## Registered S3 method overwritten by 'treeio':
## method from
## root.phylo ape
## ggtree v2.3.0.991 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use ggtree in published research, please cite the most appropriate paper(s):
##
## [36m-[39m Guangchuang Yu. Using ggtree to visualize data on tree-like structures. Current Protocols in Bioinformatics, 2020, 69:e96. doi:10.1002/cpbi.96
## [36m-[39m Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on phylogeny using ggtree. Molecular Biology and Evolution 2018, 35(12):3041-3043. doi:10.1093/molbev/msy194
## [36m-[39m Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
#读取树文件
library(treeio)
## treeio v1.11.2 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use treeio in published research, please cite:
##
## LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu. treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Molecular Biology and Evolution 2019, accepted. doi: 10.1093/molbev/msz240
#提供几何图形
library(ggstar)
#创建新的scale，多个fill或者color时
library(ggnewscale)
#为了可重复
set.seed(1024)
#读取树文件
tr <- read.tree("./tree.nwk")
#tip节点数
numtip <- length(tr$tip.label)
#dat1 有5列数据，ID列为tip名称，这一列必须映射到几何图层的y轴上, Location与Group为离散型，可以映射为几何图层的颜色或者形状
#Length与Abundance为数值型，可以映射为大小，长度等。
dat1 <- data.frame(ID=tr$tip.label,
Location=c(rep("HK", 50), rep("TW", 36), rep("SX", 30), rep("GD", 48),
rep("HN", 20), rep("AH", 20), rep("FJ", 26)),
Length=abs(rnorm(n=numtip, mean=0.6)),
Group=c(rep("Yes", 200), rep("No", 30)),
Abundance=abs(rnorm(n=numtip, mean=10, sd=3)))
dat2 <- data.frame(ID=tr$tip.label,
Pos=unlist(lapply(seq_len(8), function(x){rep(x, numtip)})),
Type=unlist(lapply(seq_len(2), function(x){rep(paste0("type",x), 4*numtip)})))
dat2 <- dat2[sample(seq_len(nrow(dat2)), 1024),]
dat3 <- data.frame(ID=rep(tr$tip.label, 4),
Type2=unlist(lapply(seq_len(4),function(x){rep(paste0("p",x), numtip)})),
Alpha=abs(rnorm(n=4*numtip, mean=0.2, sd=0.01)))
dat3$Type2 <- factor(dat3$Type2, levels=paste0("p", seq_len(4)))
dat3 <- dat3[sample(seq_len(nrow(dat3)), 200),]
# layout为树的布局，这里采用circular，size为枝的粗细，树比较大时候可以设置小点。
p <- ggtree(tr, layout="circular", size=0.1) + geom_treescale(x=6, y=0, fontsize=1.2, linesize=0.3)
# 先叠加一层树的tip节点，这里我用ggstar的图层。
# geom_fruit为ggtreeExtra提供的接口函数
p1 <- p + geom_fruit(data=dat1,
geom=geom_star,
#注意y必须映射到tip的名称列。
mapping=aes(y=ID, fill=Location, size=Length, starshape=Group),
starstroke=0.2) +
scale_size_continuous(range=c(1, 3),
guide=guide_legend(keywidth=0.5, keyheight=0.5, override.aes=list(starshape=15), order=2)) +
scale_fill_manual(values=c("#D15FEE","#EE6A50","#FFC0CB","#8E8E38","#9ACD32","#006400","#8B4513"),
guide="none")+
scale_starshape_manual(values=c(1, 15),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=1))
p1
mapping的数据可以根据自己需要展示的数据来设置。
p2 <- p1 + new_scale_fill()+
geom_fruit(data=dat2,
geom=geom_tile,
mapping=aes(y=ID, x=Pos, fill=Type),
pwidth=0.25,
position=position_identityx()) +
scale_fill_manual(values=c("#339933", "#dfac03"),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=3))
p2
fill映射在tip节点图层已经使用，所以先使用new_scale_fill()来重新初始化fill映射。而这里的position为geom_tile的参数，控制着该图层的最终位置。position_identityx函数为ggtreeExtra提供的控制整个图层的水平或者垂直平移。
p3 <- p2 + new_scale_fill() +
geom_fruit(data=dat3,
geom=geom_tile,
mapping=aes(y=ID, x=Type2, alpha=Alpha, fill=Type2),
pwidth=0.15,
position=position_identityx()) +
scale_fill_manual(values=c("#b22222", "#005500", "#0000be", "#9f1f9f"),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=4)) +
scale_alpha_continuous(range=c(0, 0.4),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=5))
p3
alpha来映射这个连续型数值。与上面一样用new_scale_fill()来重新初始化fill映射。接下来我们即将添加第四个图层–柱状图。在直角坐标系中，在树的另一端添加柱状图层很简单，ggtree已经提供了geom_facet或者facet_plot。但在极坐标的时候，这就有点难度了。其实这一点只要将整个柱状图平移下然后再转成极坐标就可以了。ggtreeExtra就是这样干的。
p4 <- p3 + new_scale_fill() +
geom_fruit(data=dat1,
geom=geom_bar,
mapping=aes(y=ID, x=Abundance, fill=Location),
pwidth=0.4,
stat="identity",
orientation="y",
position=position_stackx()) +
scale_fill_manual(values=c("#D15FEE","#EE6A50","#FFC0CB","#8E8E38","#9ACD32","#006400","#8B4513"),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=6)) +
theme(legend.position=c(0.95, 0.5),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
legend.spacing.y = unit(0.02, "cm"))
p4
legend的theme
3.2 案例2
树形图，外圈叠加柱状图与点图。上面例子外圈主要用到geom_tile与geom_bar。这个例子用到了geom_bar与geom_star。在这里我们用的树，还是上面那棵树。外圈图层的数据我们另外构造。
# Abundance为连续型，Location为离散型。
dt1 <- data.frame(ID=tr$tip.label,
Location=c(rep("HK", 50), rep("TW", 36), rep("SX", 30), rep("GD", 48),
rep("HN", 20), rep("AH", 20), rep("FJ", 26)),
Abundance=abs(rnorm(numtip, mean=20, sd=3)))
#随机选取某些tip节点的丰度设为0
dt1[sample(seq_len(nrow(dt1)), 10),"Abundance"] <- 0
dt2 <- data.frame(ID=rep(tr$tip.label,4),
Pos=rep(c(1,2), 2*numtip),
Type=rep(c("Mislabelled", "Insertions", "Refinements", "Corrections"),numtip))
dt2$Type <- factor(dt2$Type, levels=c("Mislabelled", "Insertions", "Refinements", "Corrections"))
dt2 <- dt2[sample(seq_len(nrow(dt1)), 80),]
p <- p + geom_fruit(data=dat1,
geom=geom_point,
mapping=aes(y=ID, color=Location)) +
scale_color_manual(values=c("#D15FEE","#EE6A50","#FFC0CB","#8E8E38","#9ACD32","#006400","#8B4513"),
guide="none")
f1 <- p + geom_fruit(data=dt1,
geom=geom_bar,
mapping=aes(x=Abundance, y=ID, fill=Location),
stat="identity",
orientation="y",
position=position_stackx()) +
scale_fill_manual(values=c("#D15FEE","#EE6A50","#FFC0CB","#8E8E38","#9ACD32","#006400","#8B4513"),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=1))
f1
geom_star构造的点图。
f2 <- f1 + new_scale_fill() +
geom_fruit(data=dt2,
geom=geom_star,
mapping=aes(y=ID, x=Pos, fill=Type, size=Pos),
starshape=26,
offset=0,
pwidth=0.1,
starstroke=0,
alpha=0.6,
position=position_identityx()) +
scale_fill_manual(values=c("blue", "black", "green", "red"),
guide=guide_legend(keywidth=0.5, keyheight=0.5, order=2)) +
scale_size_continuous(range=c(2,3),
guide="none") +
theme(#legend.position=c(0.96, 0.5),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
legend.spacing.y = unit(0.02, "cm"))
f2
ggstar是开发版的，才有这个starshape，由于CRAN对于版本更新频率有要求，所以还未发布的CRAN上。具体可以参考这篇《一个可以画星星图层的R包，了解下》的方法更新。
4. 总结
ggtreeExtra主要是提供了geom_fruit来将ggplot2中的图层与ggtree的tree衔接起来。具体是通过position_stackx(控制堆叠柱状图)，position_identityx(控制点，热图等图层)，position_dodgex或者position_dodgex2(控制分组柱状图，盒形图或者小提琴图)等控制图层平移，需要注意的是有些图层要用orientation参数来指定是x轴还是y轴方向。关于包名与图层函数名的确定还要感谢bibabble作者群里的小伙伴，特别要感谢李陈浩同学的建议。这个包不能画树，但是可以叠加图层到树上进行注释，类似于树上结果，所以图层就叫geom_fruit,包名就叫ggtreeExtra。
5. 往期精彩
6. 参考文献
Asnicar, Francesco, George Weingart, Timothy L Tickle, Curtis Huttenhower, and Nicola Segata. 2015. “Compact Graphical Representation of Phylogenetic Data and Metadata with Graphlan.” PeerJ 3: e1029.
Letunic, Ivica, and Peer Bork. 2019. “Interactive Tree of Life (iTOL) V4: Recent Updates and New Developments.” Nucleic Acids Research 47 (W1): W256–W259.
Subramanian, Balakrishnan, Shenghan Gao, Martin J Lercher, Songnian Hu, and Wei-Hua Chen. 2019. “Evolview V3: A Webserver for Visualization, Annotation, and Management of Phylogenetic Trees.” Nucleic Acids Research 47 (W1): W270–W275.
Yu, Guangchuang. 2020. “Using Ggtree to Visualize Data on Tree-Like Structures.” Current Protocols in Bioinformatics 69 (1): e96.https://doi.org/10.1002/cpbi.96.
Yu, Guangchuang, Tommy Tsan-Yuk Lam, Huachen Zhu, and Yi Guan. 2018. “Two Methods for Mapping and Visualizing Associated Data on Phylogeny Using Ggtree.” Molecular Biology and Evolution 35 (2): 3041–3.https://doi.org/10.1093/molbev/msy194.
Yu, Guangchuang, David Smith, Huachen Zhu, Yi Guan, and Tommy Tsan-Yuk Lam. 2017. “Ggtree: An R Package for Visualization and Annotation of Phylogenetic Trees with Their Covariates and Other Associated Data.”Methods in Ecology and Evolution 8 (1): 28–36.https://doi.org/10.1111/2041-210X.12628.